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Solving a sparse systems using linear algebra. 
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Abstract 



We generalize a method to compute the solutions of a system of polynomial equations with 
finitely many solutions. We solve this problem without the need of Grobner bases. 
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Introduction. 



In this paper we adapt a method, presented in J3, Ch. 3, §6], that solves a system with n 
polynomial equations and n variables with finitely many solutions, to a method that solves a 
general sparse system with finitely many solutions. 

In ISection II we present the general method to solve a sparse system with finitely many 
solutions. In lSection 2l and lSection 3l we adapt the method to solve a system of trilinear equations 
and a system of polynomial equations respectively. 
J> 1 The method is based on the following theorems. Consider a system of polynomial equations 

. with finitely many solutions in C", 



/i(xi,...,x„) = 



f k (xu- ■ -,x„) = 
where f\ , , . . , /j are polynomials in C[jci, . . . , x n ]- The quotient ring, 

K = C[x l ,...,x n y(f u ...,fk), 

is a finite-dimensional vector space, Theorem 2.1.2], The dimension of K is the number of 
solutions of the system counted with multiplicities. 

Every polynomial / e C[xi , . . . , x n ], determines a linear map M : 91 — > H, 

M(g)=fg, g e C[x u ...,x n ], 

where g denotes the class of the polynomial g in the quotient ring %. The matrix of M is called 
the multiplication matrix assigned to the polynomial /. 
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Theorem (Eigenvalue Theorem). The eigenvalues ofM are ),..., /(^v)}, where i, . . . , g r } 
are the solutions of the system of polynomial equations. See Theorem 2.1.4] for a proof. 

Theorem (Eigenvector Theorem). Let f = a\X\ + . . . + a n x n be a generic linear form and let M 
be its multiplication matrix. Assume that B — [l,Xi,..., X n , . . .} is a finite basis of f{ formed by 
monomials. Then the eigenvectors of M determine all the solutions of the system of polynomial 
equations. Specifically, if v — (vq, . . . ,v„, . . .) is an eigenvector of M such that vq — 1, then 
(vi, . . . , v n ) is a solution of the system of polynomial equations. Even more, every solution is of 
this form. See ^ §2.7.5/ for a proof. 

Note that the previous theorem requires that the variables {x\, . . . , x n ] are elements of the 
basis B. It could be the case that some variables are missing from B. For example, if xi, . . . , e 
B, and jc,- + i, . . . , x n B, then every missing variable, say Xj, is a linear combination of [xi, . . . , x,}, 

Xj — aj\X\ + . . . + anXi, i + 1 < j < n. 

If v = (1, Vi, V2 ■ ■ ■) is an eigenvector of M, the y'-coordinate of the solution corresponding to v, is 

0. 1 vi + . . . + ajiVj. See |6. §2.1.3]. 

The computation of the multiplication matrix is a very difficult task to do. In general, we need 
Grobner bases. When the system of polynomial equations has the same number of variables and 
equations, « = k, there exists a method to compute the multiplication matrix, see [2] and [6, 
§2.3.11. 

In yj, §7,6.21], the authors adapted the method to solve a sparse system with the same number 
of variables and equations. Here, we adapt their method to solve a general sparse system with 
finitely many solutions. We needed to prove ITheorem II in order to generalize it. The proof is 
based on standard algebro-geometric tools and Toric varieties theory. 

1. Solving a sparse system. 

A sparse system is a collection of Laurent polynomials, 

fi = £j c, '. ,,x i 1 ■ • ■ 1 - i - k > 

re<2, 

where <3, are fixed finite subsets of Z". We say that the system is generic if k = n and the 
coefficients c, ;l , are chosen arbitrarily. The set <3, is called the support of f. 

Bernstein proved in [3] that the number of solutions in (C \ {0})" of a generic sparse system 
equals the mixed volume of the corresponding Newton polytopes. The convex hull Mi of <2,, 

^ = conv(<2,0 £ K", 

is called the Newton poly tope of f, denoted N(f). 
Consider the function 

V(A U . . . , A„) := vol(AiMi + ... + A„M„), A-, > 0, 1 < i < n. 

where vol is the usual Euclidean volume in R" and M + M' denotes the Minkowski sum of 
polytopes. It is a fact that V is a homogeneous polynomial and the coefficient of the monomial 
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A\ . . .A„ is called the mixed volume of . . ., 3\ n . The mixed volume (i.e the number of solu- 
tions of a generic sparse system) is a very difficult number to compute, see [4, p. 363]. In some 

situations, this is possible and in the general case, there are a lot of algorithms to compute it. 

I — i 

When the polynomials /; are multihomogeneous, in the paper 111 111 , the author gave a recursive 
formula to compute this number. When the system is homogeneous, we recover Bezout theorem. 
Consider a sparse system with finitely many solutions, 

/l = 2 ei.vXj 1 ...x v n " = 

l'£<3i 

fk=2u Ck - vX l' ■■■ X n = 
veQk 

If the system is generic, the number of equations coincide with the number of variables, k = n, 
and the number of solutions is given by the mixed volume of the Newton poly topes of f\ , . . . , fk- 
For example, the following sparse system is generic with only one solution, £ = (1,1), 

1 - x\X2 — 

1-JC2 =0 

The supports, in Z 2 , of the equations are 

Qi = {(o,o),(i,i)}, a 2 = {(o,o),(o,i)}. 

Note that the Bezout number of this system is 2 and its mixed volume is 1 . 

Given a sparse system, it is easy to give a bound for the number of solutions. Multiply each 
equation by a positive power of a variable, and consider the corresponding Bezout number of the 
resulting system. It is always a bound for the number of solutions. If the system is a generic 
system of homogeneous polynomials, we can use Bezout's theorem to know the number of solu- 
tions, r = deg(/i) . . . deg(/„). 

Let S = C[x* 1 , . . ., x„ l ] be the algebra of Laurent polynomials. Given a Newton polytope, 
Jl, let Sjti be the vector space of polynomials with Newton polytope in M, 

= {geS : N{g) c M}. 

The dimension of S y\ is equal to the cardinal of J{ n Z", 

dim(5^) =#(^nZ"), 

The finite set J{ n Z" determines a monomial basis for Sy{. 
Let y\! be an arbitrary Newton polytope. Using the equality, 

N(gig 2 ) = N(gi) + N(g 2 ), g\,gi&S, 

we obtain that the multiplication map is surjective, 

fx : Sft®Sw — > Ssi+w, 

/Jix\ l . . . x^ <g> x} ...*?) = . . .xT v '\ Vv e M, v' e 9t. 
In other words, its image, S j\ ■ is equal to S^+^, 

Sji'Sa = S . 
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Definition. |5, Definition 2.2.9]. A Newton polytope ft is called normal if 

(k ■ ft) n z" = (ft n z") + . . . + (ft n z"), vit e n. 

The expression in the right has k terms. 

Recall from [5, Theorem 2.2. 12] that given a full dimensional Newton polytope ft, the New- 
ton polytope (n— I)- ft is normal. In particular, given ft, an arbitrary Newton polytope, (n-l)-ft 
is always normal. 

Notation. The following notations and assumptions will be used in the rest of the section. 

Let J c S be the ideal generated by a sparse system f\, . . . ,fk, where the Newton polytope 
of is ftj, 1 < i < k. 

J =</,,..., f k )Q S. 

Assume that the sparse system has finitely many solutions, say r, and that the dimension of S /ST 
is also r. In other words, all the solutions has multiplicity one. 

Let /o be a non-constant generic Laurent polynomial and consider the ideal I, 

I=(fo,fu...,fk), fo= J] c o,,^'...x);", N(f ) = fto. 

Without loss of generality, we may suppose that e ft, for every i, < i < k. Divide each 
equation by a monomial to guarantee 6 ft t , < i < k. This operation does not change the 
solutions. Recall that the solutions are in (C \ 0)". 

Assume also, that the Newton polytopes fto, . . ., ftk are normal. 

Theorem 1. Let 8 = fto + . . . + ft k and let S, •= fto + . . . + ft t + . . . + ft k , < i < k. Then the 
cokernel of the following linear map is isomorphic to S IJ~, 

k 

<t>:S Sl x...xS Sk ^S E , ®(g u ..-,gk) = Yjfi-8» 

The Newton polytope Sq C £ satisfies 

S s J(imQ>r\S So ) = Ss/im®. 
Even more, the following linear map is surjective, 

k 

Y : S So xS Bl x...xS Sl -» S 6 , x ¥(go,gu---,8k)=fo-go + Yjfi 4 Si' 

Proof. The basic references are and 

Given a very ample Newton polytope ft, J5I Definition 2.2. 17], we can construct a projective 
toric variety X and an invertible sheaf Ox(l), see J^l §8 1A]. Let X be the closure of the image 
of the following map, 

(C \ 0)" ^ P(S v ), (h,..., t n ) (?' t v »), 

where f = f"' 1 . . . C" and ft n Z" = {vi, . . . , v N ). The invertible sheaf O x (l) is the restriction of 
Op(s v )(1) to Recall from [5, Proposition 2.2. 18] that a normal Newton polytope is very ample. 
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Let Xj c P(5 „ ) be the projective toric variety assigned to the polytope J?!,-, < i < k, §5, 
IB]. Consider the following invertible sheaves in X — Xq x . . . x Xk, 

O x (d Q , ...,d k ) = n^Oxoido) ® 0x ■ ■ • ®Ox K°x k {dk), d ,...,d k eZ, 

where the map ji, ■ : X — > X/ is the /-projection, < i < k. 
By |SJ, Theorem 5.4.8], the coordinate ring of Xj is 

CK] = 0S^„ 0</<fc. 
This implies that the space of global sections of Ox(do, . . . , c4) is, 

Let us work with the projective toric variety T inside X given by, 

(?,..., f) (C \ 0)" x . . . x (C \ 0)" X = X x . . . x X k 

A i, 

t (C \ 0)" Y 

where the bottom map is the one induced by the Newton polytope £ = yio + . . . + and the 
diagram in commutative, 13, §8, 1A, Proposition 1.4]. Let Oy(do, . ..,</*) be the restriction of 
O x (do, ...,d k ) to F, 

Oy( £ /o,...,*) = O z (fl'o,...,*)|y. 

A local affine section of Ox(do, . . . ,d k ) is generated by sections of the form t^' ■ . . . ■ fT*, where 
V,- e di ■ yii and f,- = (f,i, . . . , f,„), < ; < k. Restricting to Y, we get sections f v ° + - +v *. Then 

H°(Y O Y (d , ...,d k )) = S domo+ ... +dklHk , do,..., d k > 0. 

Let e, e Z t+1 be the vector with 1 in the (z + Incoordinate and in the rest, < i < k. For 
example, eo = (1, 0, . . . , 0) and e k = (0, . . . , 0, 1). The Laurent polynomials {f\,...,fk} determine 
a <9y-linear map 

Y ^T, T = O r (ei) (£>... ®Or(ek), 
and then, we can construct the dual Koszul complex assigned to T , 

k i 

/\r w -»...-» /\r v -^...^r v ^o Y , 

where 

A^ = Oy(-e h -...-e is ). 

l<i'i<.. .<i s <k 

Let Z c (C \ 0)" be the zero locus of {/i, . . .,/*}. Given that (C \ 0)" is embedded in Y, the 
variety Z may also be embedded in Y, 

Z^Y, : ...:£•*), {vi,...,v w } = SnZ». 
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The dual Koszul complex assigned to IF, is supported in Z, §B, Proposition 1.4 (a)]. Then the 
augmented complex is exact, 

k i 

o -> /\ r v -» . . . -» /\ ^r v - . . . -» F v -» o F -» o z -» o. 

The exactness is preserved by twisting with the invertible sheaf (3y(l, . . . , 1), and taking global 
sections. In particular, the following complex of vector spaces is exact, 

S Sl X . . . X S ft A S fi -» S/J- -> 0. 

Twisting the dual Koszul complex assigned to T with the invertible sheaf Oy(0, 1, . . . , 1), we 
obtain, in a similar way, that the following map is surjective, 

S So -> S/J -» 0. 

Note that the kernel of S So -> S/J is equal to Sg n J", but given that 6 Jlo n . . . n we 
get So Q fi, and then 

S s „ n J7 = S® n S G n J7 = S So n im0> => S So /(S So n im(D) s S/J" = S 6 /im*. 

Finally, if the Laurent polynomial /o is non-zero over Z, that is, if the ideal I is equal to S , 
then similar arguments, using the sheaf T' = Oy(eo) © . . . ffi<9y(e*), show that the following map 
is surjective, 

S So x...xS St -^Ss^O. 

□ 

Consider S n V and S, n Z", < i < k, as monomial ordered bases of Ss and S s : , < i < k, 
respectively. Let p be the cardinal of So n Z", /?,• the cardinal of S, n Z", 1 < i < k and p + q be 
the cardinal of & n Z". Given that /o is non-constant, the inclusion, So £ £, is proper, thus q > 0. 

So n Z" = {mi , . . . , m p }, & n Z" = {mi , . . . , m p , m p +\, . . . , m p+9 }, 

where m, is a monomial, 1 < i < p + q. Here we are abusing the notation. The monomial 
m — Xj 1 . . . x l n " corresponds to the point (vi, . . . , v„) £ Z". 

Let M e £p+i x p+pi+--+pi< be the rectangular matrix assigned to *P in these bases, 

M = lf, U , M„ e C pxp , M 22 e C qXpi+ - +p K 

\M2l M22) 

(mi... m p m p+ \ . . . ^ j = (/ mi . . . f m p f ■ Si . . . f k ■ S k ) , 

where is the row vector obtained by multiplying /; with the monomials in S, ■ nZ", 1 < i < k. 

Theorem 2. There exists a matrix F e C Pl+ '" +PiXp such that every solution £, of the sparse 
system determine a left eigenvector of Mu + MyiF. Even more, /o(£) is the eigenvalue of that 
left eigenvector. 

The matrix F satisfies the linear condition M21 + M22F — 0. 
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Proof. First, let us see that the rank of M22 is q. The matrix M22 is the matrix of the composition 
of the following maps, 

O n 

S's, x . . . x Ss t — > Se — > Se/Ss , 
where ji is the quotient map. The rank of M22 is equal to rk(7r O), 

im(n<S>) = 7r(imO) = im<5/im(<[)) So , im(0) So = im(<$) n S So . 

Let t be the rank of M22, and let q be the dimension of SeJSg , 

t = dim(imO/im((l)) So ), q = dim(5 £ /5 So ) = dim(S fi ) - dim(S So ) = (p + q) - p. 

Using ITheorem II we obtain 

dim(5 e /imO) = dim(5 So /im(<l)) So ) => q = t. 

Now that we know that rk(M22) = q, it is easy to prove that there exists a matrix F £ C p,+ "' +PlXp 
such that 

M22F = -M 2 i. 

Each column of F, c\, . . . , c p , is a solution of the linear system M22C, = £>,, where b, e C 9 is the 
/-column vector of -M21, 1 < i < p. 

Let f e C" be a solution of the sparse system, f\(g) = . . . = = 0. Then 

(r> . . . r p r ,p+i • • • r ,p+ ") S = /o(^) • (r 1 • • • r> 0) =^ 

(f 1 . . . f *) (M, , + M 12 F) = /„<£) • (f> . . . C). 
Then, £ determines a left eigenvector of Mi 1 + M^F with eigenvalue /o(£). □ 

The following is an algorithm to compute all the solutions of a sparse system, 

2] Q.yXj 1 . . . 4" = 

Let /o be a non-constant general Laurent polynomial with N(jo) = Mo- Assume that the Newton 
polytopes jTIo, . . . , 3h are normal. In other case, take a, • Mi for some natural number a, < n - 1, 
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< ;' < k. Assume also that e yt n . . . n M k , 



Input: A sparse system with finitely many 

solutions and a non-constant Laurent polynomial. 
Output: The solutions of the system. 

L Let £ = Wo + . . . + 3Kk ™d let S, ■ = M + . . . + + . . . + 9Kk- 

2. Get ordered monomial bases, where So n Z" = {mi, . . . ,m p ] 
and £ n Z" = {mi, m p ,m p +i ... ,m p +q}. Assume m\ = 0. 

3. Compute the matrix M and decompose it as [M\\ , Mi 2; M21, M22]. 

4. Solve the linear system M22F = — M21. 

5. Compute left eigenvectors of M11+M12F. 

6. For each left eigenvector v such that vi #0do: 

6.1 Normalize v as v = (1, V2, . . .)• 

6.2 Let * be such that jc™' = v; for all 
monomial m; e 2?o n Z", 1 < ( < p. 

6.3 Save j; if it is a solution. 

7. Return the solutions. 



Remark. The formalism of Newton polytopes in our method to solve a sparse system is neces- 
sary. For example, consider the system in S = C[xi, X2], 

1 - x\X2 = 
l-x 2 =0 

Let S d be the space of polynomial of degree d. The support of a generic polynomial in S d is 

{(a l ,a 2 ) e (Z> ) 2 : a\ + a 2 < d). 

Fix a natural number e and consider the map 3), 

f :S e -2XS e -l -»5 e , 4>(5i,^2)=5i-(l-^2)+52-(l-JC2). 

Given that dim(5/(l - x x x 2 , 1 - x 2 }) = 1, we need that the codimension of im(O) be 1. But this 
never happens because the monomials 1 and x\ are linearly independent from im(<I>). 

2. Adaptation I: affine system of trilinear equations. 

Let S = C[xi, ...,x„,yi,.. .,y m ,zi, ■ ■ - ,z s ] be the polynomial algebra. Given three non- 
negatives integers, (d\ , d 2 , c/3) consider the Newton polytope 

ftd^diA = (di ■ A„) x (d 2 ■ A m ) x (d 3 ■ A. s ), 

A* = [t e (R>of : h + ... + t k <\), k &{n,m, s}. 
Let us denote S dlM instead of , 

Sd,,d 2 ,d 3 = {/ e 5 : N(f) Q Jld l ,d 2 ,d 3 }- 

We say that a polynomial in Sd u d 2 ,d 3 has multidegree less than or equal to (d\,d 2 , ^3). For exam- 
ple, the multidegree of x\x 2 y\ + ZjZ2 is less than or equal to (2, 1,3). We say that a polynomial I 
is trilinear, if the multidegree of I is less than or equal to (1, 1, 1). 
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Consider the following affine system of trilinear equations, 

' £i(x,y,z) = 

: , xeC", yeC", z e C s . 
£ k (x,y,z) = 

Assume that the system has finitely many solutions, say r. Let / be a polynomial such that, 
N(f) Q S ijj. For simplicity, take / a generic linear form. Note that /, t\, . . . , t k e S ijj. 
Bv lTheorem 11 the following linear map is surjective, 

k 

¥ : (Sk,k,kt +l — > S MM1M1 , ¥(g,gi,. ■.,8k) = g-f + Y J gi- 4 
Take B a monomial ordered basis of Sk,k,k an d extend it to a monomial ordered basis, E, of 

S k+l,k+l,k+l- 

B = 3l kM nZ" +m + s = {m u ...,m p ), 

E =^+u+u+i nZ" +m+ ' 5 = {mi m p ,m p +i, . . . ,m p+q ), 

where m,- is a monomial, 1 < i < p + q. Note that \,x\, . . . , x n ,y\, . . . ,y„„z\, ...,Zj eB. Assume 
mi = 1 ■ 

Let M e C f ' +9X(/:+1)p be the rectangular matrix in bases B and E of the linear map y ¥, 

Bv lTheorem 2l there exists a matrix F e C ipxp such that the left eigenvectors of M\\ + M\2F 

determine possible solutions of the system of trilinear equations. Even more, its eigenvalues are 

the possible values of / at the solutions of the system. 

As a corollary, we get an algorithm to compute solutions of a system of trilinear equations 
using linear algebra tools. 

Input: A system of trilinear equations with 

finitely many solutions and a generic linear form. 
Output: The solutions of the system. 

1. Construct the bases B and E. 

2. Compute the matrix M and decompose it as [M\\ , M12', M21 , M22]. 

3. Solve the linear system M22F = -Mil- 

4. Find left eigenvectors of M\\ + MiiF. 

5. For each left eigenvector v such that V] #0do: 

5.1 Normalize v as v = (1,1*2,...). 

5.2 Let x = (xi, . . . ,x n ) he such that Xj=v cr . where 
<Tj is the coordinate of x, in B, 1 < i < n. 

5.3 Same for y and z. 

5.4 Save (x, y, z) if it is a solution. 

6. Return the solutions. 



3. Adaptation II: affine system of polynomial equations. 

Consider the following affine system of polynomial equations in S = C[xu ■ ■ ■ < x n], 

/i(xi,...,x„) = 

f k (xi,...,x„) = 
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where the support of fi is {(a\, . . . , a n ) e (Z>o) n : a\+ ... + a n < di], 1 < i < k. In particular, /; 
is a polynomial of degree dj, deg(//) = dj, 1 < i < k. Assume that the system has finitely many 
solutions, say r. Let / be a generic polynomial of degree d. 

Bv lTheorem 11 there exist a degree e = d + d\ + . . . + d% such that the following linear map 
is surjective, 



*F : S e - d x S e - dl X...xS e - dk ^S e , y P(g,gu...,gk) = g-f + YjSi-fi- 



Choose monomial ordered bases of S e -d,S e -d i , . . . ,S e -a k , and extend them to a basis of S e . 
Let us denote Bq the basis of S e -d, B, the basis of S e -d t , 1 < i < k, and E the basis of S e . Let p 
be the cardinal of Bq, p, the cardinal of B,, 1 < i < k and p + q be the cardinal of E, 



where m ( - is a monomial, 1 < ; < p + q. Note that 1, x\, . . . , x n e Bo- Assume nt\ = 1. 
Let M e c p+qXp+Pl+ - +Pk be the matrix assigned to ¥, 



Bv lTheorem 21 there exists a matrix F e Cpi+- +p* x /' SU ch that the left eigenvectors of M\\ + 
MyiF determine possible solutions of the system of polynomial equations. Even more, its eigen- 
values are the possible values of / at the solutions of the system. 



Input: A system of polynomial equations with r solutions and 
a generic polynomial of degree d. 
Output: The solutions of the system. 

1. Construct the bases Bo,...,Bt,E. 

2. Compute the matrix M and decompose it as [Mi i , A/12; M21 , M22L 

3. Solve the linear system M22F = -Mil- 

4. Find left eigenvectors of M11+M12F. 

5. For each left eigenvector v such that vi # do: 

5.1 Normalize v as v = (1, v%, . . .), 

5.2 Let x = (xi, . . . ,x„) be such that jc, ■ = v tTj where 
o", is the coordinate of Xj in Bo, 1 < ' < n. 

5.3 Save x if it is a solution. 

6. Return the solutions. 



Note that Step 1 of the algorithm has exponential complexity, 0(n e ). The other steps have 
polynomial complexity. 

4. Identifying the solutions. 

Let M e (C^*^ be a diagonalizable matrix, let A e C be an eigenvalue and let v € be an 
eigenvector such that Mv = Av. In practice, there exists a positive real number s such that 



where v is the computed eigenvector corresponding to v and A is the computed eigenvalue corre- 
sponding to A. 



k 



B = {mt,...,m p }, E = {m u . . . ,m p ,m p+u . . . ,m p+q ], 




||v -v|| < e, \A- A\ < s, 
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Let x e C" be an estimated solution of the sparse system, f\ = . . . = fa = 0, obtained with 
our algorithm using the Laurent polynomial fa. Let x E C" be the exact solution corresponding 
to x~. Then, there exists a constant K depending on the coefficients of /o, . . . , /a such that, 

\fi(x)\ = \fi(x) - Mx)\ < K\\J- x\\ < K\\v- v\\ <Ks<(K+ l)e, 1 < i < k. 

\Mx) -A\< \f (x) -A\ + \A-A\< \f (x) - f (x)\ + s< KW- v|| + s < (K + l)e. 

Finally, we will identify jtas a solution of the sparse system if \ft(x)\ < (K + l)s, 1 < i < k 
and \fo(x) - A\ < (K + l)e. The constant K may be taken as the sum of the magnitudes of the 
coefficients of the polynomials. In other words, K = \\M\\\, where M is the matrix assigned to T*. 

5. Application. 

To conclude this article, let us give an application to compute the maximum of a generic 
trilinear form over a product of spheres, 

(n,m,s) 

£ : x M m+1 x R s+1 -> M, t(x,y,z)= V max \£(x,y,z)\, 

,.4^ n IWI=lly||=lk"ll=i 

(i,j,k)=0 

where the norm is the usual 2-norm. In the literature, this maximum is called the first singular 
value of {, 03 3]. 

Using Lagrange method of multipliers, |l| §13.7], the extreme points of I over a product of 
spheres, S" x x S s , satisfy 

d£/dxi(x Q ,...,x n ,y Q ,...,y m ,ZQ,...,z s ) = 2ax h < i < n, 
&e/dyj(x ,...,x n ,yo,...,ym,Zo,---,Zs) = 2/3y ; -, < j < m, 
d£/dzk(x Q ,...,x n ,y Q ,...,y m ,ZQ,...,z s ) = 2Az k , < k < s, 

a,p,AeR, \\x\\ = \\y\\ = llzll = 1. 

These equations imply that the vector d£/dx(x, y, z) is a multiple of x. Same for y and z- In other 
words, considering the system in P" x P'" x P s , we can hide the variables a,/3 and A, 

xjd£/dxi(x,y,z) = Xjd£/dxj(x,y,z), </</'< n, 
yjd£/dyi(x,y,z) = yid£/dyj(x,y,z), <i<j< m, 
Zjd£/dzi(x,y,z) = Zid£/dzj(x,y,z), <i<j< s. 

Given that £ is trilinear, the expression Xjd£/dxj(x,y,z) is equal to £(xjei,y,z), where e,- e R" +1 
is the vector with 1 in the /-coordinate and in the rest. Same for y and z- Summing up, the 
extreme points of £ satisfy the following system of equations in x R m+1 x R , 



£(xjei - Xjej,y,z) 


= o, 


< 


< ./' < n, 


£(x,yfi-yie h z) 


= o, 


< 


< /' < m 


£(x,y,Zjej -Ztej) 


= o, 


< 


< j < s, 


xl + ...+x\ 


= 1 






y + ---+y 2 m 


= 1 






zl + .-. + z] 


= 1 
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If Z is generic, we may assume that the extreme points are finite. Note that we have more equa- 
tions than variables. Enumerate the equations, /i, . . . ,fk, and let r be the number of solutions. 
We are going to apply our method to the system of polynomial equations fi,. ..,fk, and to (. Let 
X\, . . . ,A r be eigenvalues of M\\ + M\ 2 F assigned to the solutions of the system f\ — . . . — = 
and to the generic trilinear form (. If \Ai\ > 2 < i < r, then \ is the maximum value of i 
overS" x§ m xS 5 . 

The following examples were computed using the algorithm programmed in Octave, see 
|Appendix A| for a similar code. 

Example. Let f:l 2 xl 2 xl 2 ->lbe the following trilinear form 

£{x,y,z) = Ax\y\Zi + x 2 yiz\ -5xiy 2 Zi -5z\x 2 y 2 + 2x\y\z 2 -1y\x 2 z 2 -9x\y 2 z 2 -6x 2 y 2 z 2 . 
Then the maximum is, 

max \f(x,y,z)\ = 12.87128226. 

IWI=llyll=lkll=i 

Example. Let f:l 3 xl 2 xl 2 ->lbe the following trilinear form 

C{x,y,z) = 4xiyizi +x 2 yizi +6x^y\Z\ -5xiy 2 Z\ - 5 zm*2 - 6 ZW2X3+ 

2x x y\zi-l x 2 yizi - 3 x 3 y L z 2 - 9 y 2 z 2 xi - 6 y 2 z 2 x 2 - 9 y 2 z 2 x 3 . 
Then the maximum is, 

max \€(x,y,z)\ = 16.81951586. 

IMI=M=M=i 

Example. Let f:K 3 xl 3 xK 3 ->fbe the following trilinear form 

t(x,y,z) = -x 3 yiZ3 +4xiyizi +3x 3 y 2 zi -6xiy 3 z 2 -ly\Z\x 2 +4x 3 yiZi -9xiz\y 2 + 

1 Z\x 2 y 2 -5x x y 3 z\ + &x 2 y 3 zi + 2x 3 y 3 z\ + 2x\y x z 2 - 6yiz 2 x 2 - 3 x 3 y x z 2 + 
9z 2 x 2 y 2 + 3x 2 y 3 z 2 + 5 x 3 y 3 z 2 -5xiyiz 3 - 9y\z 3 x 2 - 7 xiz 3 y 2 -2x 3 y 2 z 3 +6x\y^z%+ 
7 x 2 y 3 z 3 - 10 x 3 y 3 z 3 + x x z 2 y 2 + x 3 y 2 z 2 . 

Then the maximum is, 

max \€(x,y,z)\ = 19.57534001. 

IMI=M=M=i 
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A. Octave code. 



Let us write a Octave code to solve a system of two polynomials in two variables. First of all 
we will identify a polynomial / with the following data structure, 

f(x u ...,x n )= J] c w x^...xZ\ MnZ" = {v u ...,v p }, 

f.mon{i] = Vi, f.coef{i] = c V; , 1 < i < p, /.terms = p. 
It is easy to write a program that prints a polynomial, 

function s=printPoly (f ) 
s=" " ; 

for i=l :f .terms 
v=f .mon{i} ; 
n=size(v,2) ; 
l=f .coef{i}; 
if (i — 1 1 |K0) 

s=strcat(s,num2str(l)) ; 
else 

s=strcat(s,"+",num2str(l)) ; 
end 

for j=l:n 

s=strcat(s, "*x[" ,num2str (j ) , "] "" ,num2str (v(j ) ) ) ; 
end 
end; 
end 

Another useful program is the evaluation of / at a point (xi , . . . , x„), 

function y=evalPoly (f ,x) 

y=0; 

for i=l :f .terms 

v=f .mon{i} ; 

l=f .coef{i}; 

y+=l*prod(x. ~v) ; 
end; 
end 

Finally, here is the Octave code that solves a system of two polynomials in two variables 
using fo a linear form, f\ a polynomial of degree 5 and fi of degree 7. The reader may change 
the values. 

d=[l,5,7] ; 
e=sum(d) ; 

"/Generate random polynomials of degree d_k 
for k=l : size(d, 2) 
cont=l ; 
for i=0:d(k) 
for j=0:d(k) 
if i+j<=d(k) 
f {k}.mon{cont}=[i, j] ; 
f {k} . coef {cont}=round(rand*20-10) ; 
cont++ ; 
end 
end; 
end; 

f {k} .terms=cont-l ; 
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end; 

"/■Generate monomial bases of S_{e-d_k} 
for k=l:size(d,2) 
cont=l ; 

for i=0:e-d(k) 
for j=0:e-d(k) 
if i+j<=e-d(k) 
mons{k}.pos{cont}=[i, j] ; 
cont++; 
end 
end; 
end; 

mons{k} . siz=cont-l ; 
end; 

"/Generate monomial basis of S_e appending the monomials from S_{e-d_l}. 
"/For technical reasons, the monomials are shifted by one. 
monE=zeros(e+l ,e+l) ; 
for i=l :mons{l} . siz 

v=mons{l} .pos{i}; 

monE(v(l)+l,v(2)+l)=i; 
end; 

cont=mons{l} . siz+1 ; 
for i=0:e 
for j=0:e 
if (i+j<=e && monE(i+l, j+l)==0) 
monE(i+l , j+l)=cont ; 
cont++ ; 
end 
end; 
end 

monE_siz=cont-l ; 
"/Generate Matrix M 
M=[]j 

for i=l : size (d, 2) 
for j=l:mons{i}.siz 
col=zeros(monE_siz, 1) ; 
for k=l : f{i} .terms 
v=mons{i} .pos{j}+f {i} .mon{k}; 
col(monE(v(l)+l,v(2)+l))=f{i}.coef{k}; 
end 

M=[M,col] ; 
end 
end 

"/Generate Matrix decomposition, 
M11=M(1 :mons{l} . siz, 1 :mons{l} . siz) ; 
M12=M(1 :mons{l} . siz,mons{l} . siz+1 : size (M, 2) ) ; 
M21=M(mons{l} . siz+1 : size(M, 1) , 1 :mons{l} . siz) ; 
M22=M(mons{l} . siz+1 : size (M, 1) ,mons{l} . siz+1 : size(M, 2) ) ; 
"/Compute left eigenvectors 
F=M22\(-M21) ; 
res=Mll+M12*F; 
[v,l] =eig(res ' ) ; 
'/Find solutions 
for k=l : size (v,2) 
if abs(v(l,k))>0 

posxl=monE(2, 1) ; 

posx2=monE(l,2) ; 

x= [v(posxl,k) ,v(posx2,k)]/v(l,k) ; 

if norm([evalPoly(f{l},x)-l(k,k) , evalPoly (f {2} ,x) ,evalPoly(f{3},x)])<l 
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X 

end 
end 
end 

The following table shows the time in seconds used to solve Step 1 throw 4 of the algorithm. 
In the first column appears the degree of the polynomials, in the second the time of our method 
and in the third, the time used with a Grobner Basis algorithm. All the computations were made 
on a 2.1GHz CPU, with 2GB of memory. We used part of the previous code to generate the 
second column and we used the following Maple 1 1 code to generate the third column, 

> G:=Basis([f [2] ,f [3]] , 'tord') : 

> ns,rv:=NormalSet(G, tord): 

> mulMat : =MultiplicationMatrix(f [1] , ns ,rv,G ,tord) : 

> evalf (eigenvectors (evalf (mulMat) )) : 

The resulting table is the following, 



degrees 


Steps 1-4 


Grobner 


(1,5,7) 


0.24 


2.56 


(1,5,9) 


0.37 


4.91 


(1,5,11) 


0.56 


9.77 


(1,7,9) 


0.62 


14.46 


(1,7,11) 


0.89 


25.19 


(1,9,11) 


1.34 


57.67 


(1,11,11) 


2.04 


114.24 
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